function [FSL]=codeFSL
 load mydata;
% cd('E:');
N=length(X1);
M=19;  %******%
%% refrence state & W1, W2
% central point
if N-2*fix((N/2))==0
   Nref=N/2;
else
    Nref=(N+1)/2;
end
% refrence state for each channel
Xref(1,:)=X1(Nref,:);
Xref(2,:)=X2(Nref,:);
Xref(3,:)=X3(Nref,:);
Xref(4,:)=X4(Nref,:);
Xref(5,:)=X5(Nref,:); 
Xref(6,:)=X6(Nref,:); 
Xref(7,:)=X7(Nref,:);
Xref(8,:)=X8(Nref,:);
Xref(9,:)=X9(Nref,:); 
Xref(10,:)=X10(Nref,:);
Xref(11,:)=X11(Nref,:);
Xref(12,:)=X12(Nref,:); 
Xref(13,:)=X13(Nref,:); 
Xref(14,:)=X14(Nref,:);
Xref(15,:)=X15(Nref,:);
Xref(16,:)=X16(Nref,:);
Xref(17,:)=X17(Nref,:);
Xref(18,:)=X18(Nref,:);
Xref(19,:)=X19(Nref,:);
 % deletion refrence state from each channel   
X1(Nref,:)=[];
X2(Nref,:)=[];
X3(Nref,:)=[];
X4(Nref,:)=[];
X5(Nref,:)=[];   
X6(Nref,:)=[];
X7(Nref,:)=[];
X8(Nref,:)=[];
X9(Nref,:)=[];
X10(Nref,:)=[];
X11(Nref,:)=[];
X12(Nref,:)=[];
X13(Nref,:)=[];    
X14(Nref,:)=[];    
X15(Nref,:)=[];    
X16(Nref,:)=[];    
X17(Nref,:)=[];    
X18(Nref,:)=[];    
X19(Nref,:)=[];    
N=length(X1);     
% Euclidean distance between states & refrence state
for i=1:N
 dis1(i)=Euclidean(X1(i,:),Xref(1,:));
 dis2(i)=Euclidean(X2(i,:),Xref(2,:));
 dis3(i)=Euclidean(X3(i,:),Xref(3,:));
 dis4(i)=Euclidean(X4(i,:),Xref(4,:));
 dis5(i)=Euclidean(X5(i,:),Xref(5,:));
 dis6(i)=Euclidean(X6(i,:),Xref(6,:));
 dis7(i)=Euclidean(X7(i,:),Xref(7,:));
 dis8(i)=Euclidean(X8(i,:),Xref(8,:));
 dis9(i)=Euclidean(X9(i,:),Xref(9,:));
 dis10(i)=Euclidean(X10(i,:),Xref(10,:));
 dis11(i)=Euclidean(X11(i,:),Xref(11,:));
 dis12(i)=Euclidean(X12(i,:),Xref(12,:));
 dis13(i)=Euclidean(X13(i,:),Xref(13,:));
 dis14(i)=Euclidean(X14(i,:),Xref(14,:));
 dis15(i)=Euclidean(X15(i,:),Xref(15,:));
 dis16(i)=Euclidean(X16(i,:),Xref(16,:));
 dis17(i)=Euclidean(X17(i,:),Xref(17,:));
 dis18(i)=Euclidean(X18(i,:),Xref(18,:));
 dis19(i)=Euclidean(X19(i,:),Xref(19,:));
end
% compute w1 & w2, w1 is min distance of a state to refrence
% and w2 is max distance
w1(1)=min(dis1);  w1(2)=min(dis2);  w1(3)=min(dis3);  w1(4)=min(dis4);    w1(5)=min(dis5);    w1(6)=min(dis6);
w1(7)=min(dis7);  w1(8)=min(dis8);  w1(9)=min(dis9);  w1(10)=min(dis10);  w1(11)=min(dis11);  w1(12)=min(dis12);
w1(13)=min(dis13);w1(14)=min(dis14);w1(15)=min(dis15);w1(16)=min(dis16);w1(17)=min(dis17);w1(18)=min(dis18);w1(19)=min(dis19);
w2(1)=max(dis1);  w2(2)=max(dis2);  w2(3)=max(dis3);  w2(4)=max(dis4);    w2(5)=max(dis5);    w2(6)=max(dis6);  
w2(7)=max(dis7);  w2(8)=max(dis8);  w2(9)=max(dis9);  w2(10)=max(dis10);  w2(11)=max(dis11);  w2(12)=max(dis12);
w2(13)=max(dis13); w2(14)=max(dis14);w2(15)=max(dis15);w2(16)=max(dis16);w2(17)=max(dis17);w2(18)=max(dis18);w2(19)=max(dis19);
%% optimum epsilon
% we must choose epsilon in a way that the probability of nearing 
% of trajectory points equals Pref.
%eps=[163,454,128,257,812,136,304,563,124,274,424,133];
eps0=0.1*ones(1,19);   % initial epsilon
Pref=0.05;   % refrence P
eps=optimum_epsilon_FSL(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17,X18,X19,Xref,w1,w2,eps0,Pref)
%% formol 1
for m=1:N
mu(1,m)=exp(-(Euclidean(Xref(1,:),X1(m,:))^2)/eps(1));
end
for m=1:N
mu(2,m)=exp(-(Euclidean(Xref(2,:),X2(m,:))^2)/eps(2));
end
for m=1:N
mu(3,m)=exp(-(Euclidean(Xref(3,:),X3(m,:))^2)/eps(3));
end
for m=1:N
mu(4,m)=exp(-(Euclidean(Xref(4,:),X4(m,:))^2)/eps(4));
end
for m=1:N
mu(5,m)=exp(-(Euclidean(Xref(5,:),X5(m,:))^2)/eps(5));
end
for m=1:N
mu(6,m)=exp(-(Euclidean(Xref(6,:),X6(m,:))^2)/eps(6));
end
for m=1:N
mu(7,m)=exp(-(Euclidean(Xref(7,:),X7(m,:))^2)/eps(7));
end
for m=1:N
mu(8,m)=exp(-(Euclidean(Xref(8,:),X8(m,:))^2)/eps(8));
end
for m=1:N
mu(9,m)=exp(-(Euclidean(Xref(9,:),X9(m,:))^2)/eps(9));
end
for m=1:N
mu(10,m)=exp(-(Euclidean(Xref(10,:),X10(m,:))^2)/eps(10));
end
for m=1:N
mu(11,m)=exp(-(Euclidean(Xref(11,:),X11(m,:))^2)/eps(11));
end
for m=1:N
mu(12,m)=exp(-(Euclidean(Xref(12,:),X12(m,:))^2)/eps(12));
end
for m=1:N
mu(13,m)=exp(-(Euclidean(Xref(13,:),X13(m,:))^2)/eps(13));
end
for m=1:N
mu(14,m)=exp(-(Euclidean(Xref(14,:),X14(m,:))^2)/eps(14));
end
for m=1:N
mu(15,m)=exp(-(Euclidean(Xref(15,:),X15(m,:))^2)/eps(15));
end
for m=1:N
mu(16,m)=exp(-(Euclidean(Xref(16,:),X16(m,:))^2)/eps(16));
end
for m=1:N
mu(17,m)=exp(-(Euclidean(Xref(17,:),X17(m,:))^2)/eps(17));
end
for m=1:N
mu(18,m)=exp(-(Euclidean(Xref(18,:),X18(m,:))^2)/eps(18));
end
for m=1:N
mu(19,m)=exp(-(Euclidean(Xref(19,:),X19(m,:))^2)/eps(19));
end
 %% formol 3
 H(m)=0;
 for m=1:N
 H(m)=sum(mu(1:end,m))+H(m);%sum(mu(1:end,m));
 end                                  
 %% formol 4
for k=1:M
     W(k)=0;
 for m=1:N
     W(k)=W(k)+(((H(m)-mu(k,m))/(M-1))*mu(k,m));
 end
end

for k=1:M
   sum2(k)=0;
   for m=1:N
    sum2(k)=sum2(k)+mu(k,m);
   end
end

% fuzzy synchronization likelihood of 19 channels 
for k=1:M
    FSL(k)=W(k)./sum2(k);
end
 
 